Regression (Phase-based) CLEAN

Intro

Here we provide exploration of LMER models fitted to the 2020 U.S. Congressional election Twitter Dataset (https://doi.org/10.6084/m9.figshare.25523734).

This document was typeset assembled using Quarto (https://quarto.org).

Initialize Workspace

Clear and Init Workspace

Imports

library(ggplot2)
library(GGally)
library(plyr)
library(tidyverse)

Set Directories

Data Load and Inspection

Data Load

ifn0 <- "i0001-politicians-aggregated.rds"
suppressWarnings(rm(list = ls(pattern = "^df")))
df0 <- readr::read_rds(file.path(ifd0, ifn0)) %>%
  dplyr::select(Name, Party, Outcome, Period, Phase, Time, Days, Agency) %>%
  dplyr::mutate(
    Period = forcats::fct_recode(
      Period,
      PC = "MI",
      PB = "BE",
      PA = "AR")) %>%
  identity()
df0 %>% dim() %>% cat0()

169997 8 

for (name in names(df0)) {cat0(paste0("    ", name, ", "))}

    Name,  
    Party,  
    Outcome,  
    Period,  
    Phase,  
    Time,  
    Days,  
    Agency,  

Skim

set.w(222)
library(skimr)
options(scipen=999)
skimr::skim_format(numeric=list(digits=6))

Warning: Formatting options are now in print() or as function arguments in skim_with().

## digits=4, scientific = FALSE
df0 %>% skimr::skim() %>% dplyr::select(-c(n_missing)) %>% print(numeric=list(digits=2))

── Data Summary ────────────────────────
                           Values    
Name                       Piped data
Number of rows             169997    
Number of columns          8         
_______________________              
Column type frequency:               
  factor                   5         
  numeric                  3         
________________________             
Group variables            None      

── Variable type: factor ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable complete_rate ordered n_unique top_counts                                
1 Name                      1 FALSE        870 Bar: 361, Ela: 361, Llo: 361, Pra: 361    
2 Party                     1 FALSE          2 Dem: 99081, Rep: 70916                    
3 Outcome                   1 FALSE          2 win: 115031, los: 54966                   
4 Period                    1 FALSE          3 PB: 98660, PA: 47974, PC: 23363           
5 Phase                     1 FALSE          4 BE: 98660, AR: 47974, AE: 12011, BR: 11352

── Variable type: numeric ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable complete_rate     mean      sd      p0     p25     p50    p75   p100 hist 
1 Time                      1  -0.0680   0.572   -1     -0.55   -0.122  0.428   1    ▇▇▇▆▆
2 Days                      1 -12.2    103.    -180    -99     -22     77     180    ▇▇▇▆▆
3 Agency                    1   0.500    0.271   -1.69   0.343   0.500  0.658   2.37 ▁▁▇▂▁

Check factors order

df2 <- df0
df3 <- df0
summary(df0$Period)

   PC    PB    PA 
23363 98660 47974 

df2$Period <- factor(df2$Period, levels = c("PB", "PC", "PA"))
summary(df2$Period)

   PB    PC    PA 
98660 23363 47974 

df3$Period <- relevel(df0$Period, ref = "PB")
summary(df3$Period)

   PB    PC    PA 
98660 23363 47974 

Phase Levels

summary(df2$Phase) ##  <- factor(df2$Phase, levels = c("PB", "PC", "PA"))

   BE    AE    BR    AR 
98660 12011 11352 47974 

attributes(df2$Phase)

$levels
[1] "BE" "AE" "BR" "AR"

$class
[1] "factor"

Check Phase and Period

df0 %>%
  dplyr::group_by(Phase, Period) %>%
  dplyr::summarize(
    DaysMin = min(Days),
    DaysMax = max(Days),
    .groups = "drop")

# A tibble: 4 × 4
  Phase Period DaysMin DaysMax
  <fct> <fct>    <dbl>   <dbl>
1 BE    PB        -180       0
2 AE    PC           1      32
3 BR    PC          33      63
4 AR    PA          64     180

Prepare for Regression

File names and output directories

## ./data/j1002-model-by-period/i0001_models/fit0x0.extension
## ./data/j1008-model-by-phase/i0001_models/fit0x0.extension

Stuff for plotting

Marginalization

Checkups

Save data frame for reference

file <- file.path(dir8, "data-df0.rds")
df0 %>% readr::write_rds(file)
cat0(file)

./data/j1008-model-by-phase/i0001_models/data-df0.rds 

Model 01: Null

Fit

model <- "fit01aPh"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (1 | Name) + 1,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)

fit01aPh: [df0] Agency ~ (1 | Name) + 1

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (1 | Name) + 1
   Data: df0
Control: control

REML criterion at convergence: 26631.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-7.9514 -0.5639 -0.0023  0.5683  7.3986 

Random effects:
 Groups   Name        Variance Std.Dev.
 Name     (Intercept) 0.006209 0.0788  
 Residual             0.067519 0.2598  
Number of obs: 169997, groups:  Name, 870

Fixed effects:
              Estimate Std. Error         df t value            Pr(>|t|)    
(Intercept)   0.498618   0.002807 826.426886   177.7 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat0(sep2); pp(model); performance::r2(get(model))

--------------------------------------------------------------------- 

fit01aPh: [df0] Agency ~ (1 | Name) + 1

# R2 for Mixed Models

  Conditional R2: 0.084
     Marginal R2: 0.000

cat0(sep2); pp(model); performance::icc(get(model))

--------------------------------------------------------------------- 

fit01aPh: [df0] Agency ~ (1 | Name) + 1

# Intraclass Correlation Coefficient

    Adjusted ICC: 0.084
  Unadjusted ICC: 0.084

cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)

--------------------------------------------------------------------- 

fit01aPh: [df0] Agency ~ (1 | Name) + 1

# ICC by Group

Group |   ICC
-------------
Name  | 0.084

Random Effects

extra <- "9001"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="re")

file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
gg88

Model 02: Time

Fit

model <- "fit02aPh"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (Time | Name) + Time,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)

fit02aPh: [df0] Agency ~ (Time | Name) + Time

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (Time | Name) + Time
   Data: df0
Control: control

REML criterion at convergence: 24941.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.0353 -0.5639 -0.0047  0.5664  7.3362 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Name     (Intercept) 0.006283 0.07927      
          Time        0.004088 0.06394  0.18
 Residual             0.066384 0.25765      
Number of obs: 169997, groups:  Name, 870

Fixed effects:
              Estimate Std. Error         df t value             Pr(>|t|)    
(Intercept)   0.498029   0.002848 816.443563 174.884 < 0.0000000000000002 ***
Time         -0.010230   0.002634 709.116999  -3.884             0.000112 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
     (Intr)
Time 0.187 

cat0(sep2); pp(model); performance::r2(get(model))

--------------------------------------------------------------------- 

fit02aPh: [df0] Agency ~ (Time | Name) + Time

# R2 for Mixed Models

  Conditional R2: 0.102
     Marginal R2: 0.000

cat0(sep2); pp(model); performance::icc(get(model))

--------------------------------------------------------------------- 

fit02aPh: [df0] Agency ~ (Time | Name) + Time

# Intraclass Correlation Coefficient

    Adjusted ICC: 0.102
  Unadjusted ICC: 0.102

cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)

--------------------------------------------------------------------- 

fit02aPh: [df0] Agency ~ (Time | Name) + Time

Model contains random slopes. Cannot compute accurate ICCs by group factors.

# ICC by Group

Group |   ICC
-------------
Name  | 0.085

Model Summary

## pp(model);
tbl0 <- gtsummary::tbl_regression(
  get(model),
  add_pairwise_contrasts = TRUE,
  exponentiate = FALSE,
  tidy_fun = broom.mixed::tidy,
)
tbl0

Characteristic	Beta	95% CI1	p-value
Time	-0.01	-0.02, -0.01	<0.001
Name.sd__(Intercept)	0.08
Name.cor__(Intercept).Time	0.18
Name.sd__Time	0.06
Residual.sd__Observation	0.26
1 CI = Confidence Interval

Model Parameters

pp(model);

fit02aPh: [df0] Agency ~ (Time | Name) + Time

parameters::model_parameters(get(model))

# Fixed Effects

Parameter   | Coefficient |       SE |         95% CI | t(169991) |      p
--------------------------------------------------------------------------
(Intercept) |        0.50 | 2.85e-03 | [ 0.49,  0.50] |    174.88 | < .001
Time        |       -0.01 | 2.63e-03 | [-0.02, -0.01] |     -3.88 | < .001

# Random Effects

Parameter                  | Coefficient
----------------------------------------
SD (Intercept: Name)       |        0.08
SD (Time: Name)            |        0.06
Cor (Intercept~Time: Name) |        0.18
SD (Residual)              |        0.26

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.

Plot Model: Estimates

extra <- "0000"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="est") 
gg88 <- gg88 + ggplot2::ylim(-0.1,0.1) + line2 + stuff0

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time

terms <- c("Time")
extra <- "0001"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

EF: Time

NULL: Time

extra <- "1001"
terms <- c("Time")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

pp(model);

fit02aPh: [df0] Agency ~ (Time | Name) + Time

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.51 | 0.50, 0.51
-0.50 |      0.50 | 0.50, 0.51
 0.00 |      0.50 | 0.49, 0.50
 0.50 |      0.50 | 0.49, 0.50
 1.00 |      0.49 | 0.48, 0.50

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Slope     |       95% CI |     p
--------------------------------
-7.57e-03 | -0.01,  0.00 | 0.004

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.15,)

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0
## gg88 <- gg88 + coord_cartesian(ylim = c(0.45, 0.55))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

Random Effects

extra <- "9001"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="re")

file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
gg88

Model 03: Time x Phase

Fit

model <- "fit03aPh"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (Time | Name) + Time * Phase,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (Time | Name) + Time * Phase
   Data: df0
Control: control

REML criterion at convergence: 24187.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.1124 -0.5632 -0.0069  0.5646  7.3739 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Name     (Intercept) 0.006236 0.07897      
          Time        0.004124 0.06422  0.18
 Residual             0.066073 0.25705      
Number of obs: 169997, groups:  Name, 870

Fixed effects:
                  Estimate    Std. Error            df t value             Pr(>|t|)    
(Intercept)       0.528553      0.003172   1283.400354 166.622 < 0.0000000000000002 ***
Time              0.044610      0.003657   2775.135404  12.198 < 0.0000000000000002 ***
PhaseAE           0.015184      0.004861 168870.559183   3.123              0.00179 ** 
PhaseBR          -0.132158      0.012989 168617.283756 -10.175 < 0.0000000000000002 ***
PhaseAR          -0.000602      0.004747 169706.699223  -0.127              0.89909    
Time:PhaseAE     -0.638482      0.044797 168682.837032 -14.253 < 0.0000000000000002 ***
Time:PhaseBR      0.224563      0.047871 168575.449036   4.691           0.00000272 ***
Time:PhaseAR     -0.110523      0.007010 167010.782296 -15.765 < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) Time   PhasAE PhasBR PhasAR Tm:PAE Tm:PBR
Time         0.431                                          
PhaseAE     -0.165 -0.219                                   
PhaseBR     -0.061 -0.083  0.043                            
PhaseAR     -0.171 -0.223  0.117  0.044                     
Time:PhasAE -0.026 -0.049 -0.786  0.009  0.022              
Time:PhasBR -0.026 -0.045  0.017 -0.966  0.019  0.004       
Time:PhasAR -0.172 -0.314  0.117  0.046 -0.695  0.030  0.025

cat0(sep2); pp(model); performance::r2(get(model))

--------------------------------------------------------------------- 

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

# R2 for Mixed Models

  Conditional R2: 0.106
     Marginal R2: 0.005

cat0(sep2); pp(model); performance::icc(get(model))

--------------------------------------------------------------------- 

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

# Intraclass Correlation Coefficient

    Adjusted ICC: 0.102
  Unadjusted ICC: 0.101

cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)

--------------------------------------------------------------------- 

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

Model contains random slopes. Cannot compute accurate ICCs by group factors.

# ICC by Group

Group |   ICC
-------------
Name  | 0.085

Model Summary

pp(model);

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

tbl0 <- gtsummary::tbl_regression(
  get(model),
  add_pairwise_contrasts = TRUE,
  exponentiate = FALSE,
  tidy_fun = broom.mixed::tidy,
)

Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 169997' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.

Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 169997' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.

NOTE: Results may be misleading due to involvement in interactions

tbl0

Characteristic	Beta	95% CI1	p-value
Time	0.04	0.04, 0.05	<0.001
Phase
AE - BE	0.06	0.04, 0.08	<0.001
BR - BE	-0.15	-0.19, -0.11	<0.001
BR - AE	-0.21	-0.25, -0.16	<0.001
AR - BE	0.01	-0.01, 0.02	0.5
AR - AE	-0.05	-0.07, -0.03	<0.001
AR - BR	0.15	0.11, 0.20	<0.001
Time * Phase
Time * AE	-0.64	-0.73, -0.55	<0.001
Time * BR	0.22	0.13, 0.32	<0.001
Time * AR	-0.11	-0.12, -0.10	<0.001
Name.sd__(Intercept)	0.08
Name.cor__(Intercept).Time	0.18
Name.sd__Time	0.06
Residual.sd__Observation	0.26
1 CI = Confidence Interval

Model Parameters

pp(model);

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

parameters::model_parameters(get(model))

# Fixed Effects

Parameter         | Coefficient |       SE |         95% CI | t(169985) |      p
--------------------------------------------------------------------------------
(Intercept)       |        0.53 | 3.17e-03 | [ 0.52,  0.53] |    166.62 | < .001
Time              |        0.04 | 3.66e-03 | [ 0.04,  0.05] |     12.20 | < .001
Phase [AE]        |        0.02 | 4.86e-03 | [ 0.01,  0.02] |      3.12 | 0.002 
Phase [BR]        |       -0.13 |     0.01 | [-0.16, -0.11] |    -10.17 | < .001
Phase [AR]        |   -6.02e-04 | 4.75e-03 | [-0.01,  0.01] |     -0.13 | 0.899 
Time × Phase [AE] |       -0.64 |     0.04 | [-0.73, -0.55] |    -14.25 | < .001
Time × Phase [BR] |        0.22 |     0.05 | [ 0.13,  0.32] |      4.69 | < .001
Time × Phase [AR] |       -0.11 | 7.01e-03 | [-0.12, -0.10] |    -15.77 | < .001

# Random Effects

Parameter                  | Coefficient
----------------------------------------
SD (Intercept: Name)       |        0.08
SD (Time: Name)            |        0.06
Cor (Intercept~Time: Name) |        0.18
SD (Residual)              |        0.26

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.

Plot Model: Estimates

pp(model)

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

extra <- "0000"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="est")
gg88 <- gg88 + ggplot2::ylim(-0.3,0.3) + line2 + stuff0

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).

file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).

gg88

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).

FE: Time

pp(model);

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

extra <- "0001"
terms=c("Time")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Phase

pp(model);

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

extra <- "0002"
terms=c("Phase")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0


file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Phase

pp(model);

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

extra <- "0003"
terms=c("Time", "Phase")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + lineR + line1 + stuff0 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

EF: Time

NULL: Time

pp(model)

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

extra <- "1001"
terms <- c("Time")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.54 | 0.52, 0.55
-0.50 |      0.53 | 0.52, 0.54
 0.00 |      0.52 | 0.52, 0.53
 0.50 |      0.52 | 0.51, 0.52
 1.00 |      0.51 | 0.50, 0.52

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Slope |       95% CI |     p
----------------------------
-0.01 | -0.02,  0.00 | 0.025

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05,) 
gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 
## gg88 <- gg88 + coord_cartesian(ylim = c(0.25, 0.75))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

EF: Phase

NULL: Phase

pp(model);

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

extra <- "1002"
terms <- c("Phase")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Phase | Predicted |     95% CI
------------------------------
BE    |      0.53 | 0.52, 0.53
AE    |      0.59 | 0.57, 0.60
BR    |      0.38 | 0.35, 0.41
AR    |      0.54 | 0.52, 0.55

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

Phase | Predicted |     95% CI |      p
---------------------------------------
BE    |      0.53 | 0.52, 0.53 | < .001
AE    |      0.59 | 0.57, 0.60 | < .001
BR    |      0.38 | 0.35, 0.41 | < .001
AR    |      0.54 | 0.52, 0.55 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Phase

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
## gg88 <- gg88 + line0 + line1 + time0 + stuff0 

gg88 <- gg88 + stuff0


file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Phase

pp(model);

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# Pairwise comparisons

Phase |  Contrast |       95% CI |      p
-----------------------------------------
BE-AE |     -0.06 | -0.07, -0.04 | < .001
BE-BR |      0.15 |  0.12,  0.18 | < .001
BE-AR | -6.92e-03 | -0.02,  0.00 | 0.174 
AE-BR |      0.21 |  0.17,  0.24 | < .001
AE-AR |      0.05 |  0.03,  0.07 | < .001
BR-AR |     -0.15 | -0.19, -0.12 | < .001

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Phase

NULL: Time x Phase

pp(model);

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

extra <- "1003"
terms <- c("Time", "Phase")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")


cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Phase: BE

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.48 | 0.47, 0.49
-0.50 |      0.51 | 0.50, 0.51
 0.00 |      0.53 | 0.52, 0.54
 0.50 |      0.56 | 0.55, 0.56
 1.00 |      0.58 | 0.57, 0.59

Phase: AE

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      1.14 |  1.04, 1.23
-0.50 |      0.84 |  0.79, 0.89
 0.00 |      0.55 |  0.54, 0.56
 0.50 |      0.25 |  0.21, 0.29
 1.00 |     -0.04 | -0.12, 0.04

Phase: BR

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.13 | 0.01, 0.24
-0.50 |      0.26 | 0.19, 0.33
 0.00 |      0.40 | 0.37, 0.42
 0.50 |      0.54 | 0.51, 0.56
 1.00 |      0.67 | 0.60, 0.74

Phase: AR

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.59 | 0.57, 0.61
-0.50 |      0.56 | 0.55, 0.58
 0.00 |      0.53 | 0.52, 0.54
 0.50 |      0.50 | 0.49, 0.51
 1.00 |      0.47 | 0.46, 0.48

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Phase | Slope |       95% CI |      p
-------------------------------------
BE    |  0.04 |  0.04,  0.05 | < .001
AE    | -0.59 | -0.68, -0.51 | < .001
BR    |  0.27 |  0.18,  0.36 | < .001
AR    | -0.07 | -0.08, -0.05 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time x Period

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 
gg88 <- gg88 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Time x Phase

pp(model);

fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# (Average) Linear trend for Time

Phase | Contrast |       95% CI |      p
----------------------------------------
BE-AE |     0.64 |  0.55,  0.73 | < .001
BE-BR |    -0.22 | -0.32, -0.13 | < .001
BE-AR |     0.11 |  0.10,  0.12 | < .001
AE-BR |    -0.86 | -0.99, -0.73 | < .001
AE-AR |    -0.53 | -0.62, -0.44 | < .001
BR-AR |     0.34 |  0.24,  0.43 | < .001

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

Model 04: Time x Phase x Outcome

Fit

model <- "fit04aPh"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (Time | Name) + Time * Phase * Outcome,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (Time | Name) + Time * Phase * Outcome
   Data: df0
Control: control

REML criterion at convergence: 23685.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.1858 -0.5634 -0.0081  0.5638  7.3838 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 Name     (Intercept) 0.004804 0.06931       
          Time        0.003343 0.05782  -0.06
 Residual             0.065972 0.25685       
Number of obs: 169997, groups:  Name, 870

Fixed effects:
                                Estimate    Std. Error            df t value             Pr(>|t|)    
(Intercept)                     0.515092      0.004452   1470.335766 115.699 < 0.0000000000000002 ***
Time                            0.057026      0.005598   3432.607484  10.188 < 0.0000000000000002 ***
PhaseAE                        -0.017658      0.008833 169175.854063  -1.999              0.04559 *  
PhaseBR                        -0.145312      0.026930 168955.292605  -5.396         0.0000000683 ***
PhaseAR                        -0.104718      0.009831 166212.881574 -10.652 < 0.0000000000000002 ***
Outcomewinner                   0.023922      0.005856   1424.761477   4.085         0.0000465435 ***
Time:PhaseAE                   -0.882126      0.086357 169124.880897 -10.215 < 0.0000000000000002 ***
Time:PhaseBR                    0.017444      0.099141 168923.961312   0.176              0.86033    
Time:PhaseAR                   -0.077030      0.014696 150506.476445  -5.241         0.0000001596 ***
Time:Outcomewinner             -0.020359      0.007198   3282.626151  -2.828              0.00471 ** 
PhaseAE:Outcomewinner           0.051078      0.010583 169008.316334   4.826         0.0000013927 ***
PhaseBR:Outcomewinner           0.022108      0.030740 168862.263648   0.719              0.47202    
PhaseAR:Outcomewinner           0.140686      0.011241 167976.986518  12.516 < 0.0000000000000002 ***
Time:PhaseAE:Outcomewinner      0.306640      0.101049 168997.306098   3.035              0.00241 ** 
Time:PhaseBR:Outcomewinner      0.281534      0.113190 168837.528411   2.487              0.01287 *  
Time:PhaseAR:Outcomewinner     -0.035756      0.016747 157881.162351  -2.135              0.03276 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of fixed effects could have been required in summary()

Correlation of Fixed Effects:
            (Intr) Time   PhasAE PhasBR PhasAR Otcmwn Tm:PAE Tm:PBR Tm:PAR Tm:Otc PhAE:O PhBR:O PhAR:O T:PAE: T:PBR:
Time         0.364                                                                                                  
PhaseAE     -0.156 -0.194                                                                                           
PhaseBR     -0.051 -0.064  0.029                                                                                    
PhaseAR     -0.141 -0.177  0.079  0.026                                                                             
Outcomewnnr -0.760 -0.277  0.119  0.039  0.108                                                                      
Time:PhasAE -0.024 -0.043 -0.783  0.008  0.017  0.018                                                               
Time:PhasBR -0.021 -0.037  0.011 -0.973  0.013  0.016  0.002                                                        
Time:PhasAR -0.145 -0.248  0.075  0.027 -0.794  0.110  0.023  0.016                                                 
Tm:Otcmwnnr -0.283 -0.778  0.151  0.050  0.138  0.355  0.033  0.029  0.193                                          
PhsAE:Otcmw  0.130  0.162 -0.835 -0.024 -0.066 -0.169  0.653 -0.009 -0.063 -0.211                                   
PhsBR:Otcmw  0.045  0.056 -0.025 -0.876 -0.023 -0.058 -0.007  0.852 -0.024 -0.073  0.034                            
PhsAR:Otcmw  0.124  0.155 -0.069 -0.023 -0.875 -0.161 -0.015 -0.012  0.694 -0.200  0.094  0.032                     
Tm:PhsAE:Ot  0.020  0.037  0.669 -0.007 -0.014 -0.026 -0.855 -0.002 -0.020 -0.047 -0.783  0.008  0.018              
Tm:PhsBR:Ot  0.019  0.032 -0.010  0.852 -0.012 -0.024 -0.002 -0.876 -0.014 -0.041  0.014 -0.970  0.015  0.003       
Tm:PhsAR:Ot  0.127  0.218 -0.066 -0.024  0.696 -0.162 -0.021 -0.014 -0.878 -0.280  0.092  0.034 -0.759  0.025  0.019

cat0(sep2); pp(model); performance::r2(get(model))

--------------------------------------------------------------------- 

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

# R2 for Mixed Models

  Conditional R2: 0.101
     Marginal R2: 0.020

cat0(sep2); pp(model); performance::icc(get(model))

--------------------------------------------------------------------- 

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

# Intraclass Correlation Coefficient

    Adjusted ICC: 0.083
  Unadjusted ICC: 0.081

cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)

--------------------------------------------------------------------- 

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

Model contains random slopes. Cannot compute accurate ICCs by group factors.

# ICC by Group

Group |   ICC
-------------
Name  | 0.067

Model Summary

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

tbl0 <- gtsummary::tbl_regression(
  get(model),
  add_pairwise_contrasts = TRUE,
  exponentiate = FALSE,
  tidy_fun = broom.mixed::tidy,
)

Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 169997' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.

Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 169997' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.

NOTE: Results may be misleading due to involvement in interactions

Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 169997' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.

Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 169997' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.

NOTE: Results may be misleading due to involvement in interactions

tbl0

Characteristic	Beta	95% CI1	p-value
Time	0.06	0.05, 0.07	<0.001
Phase
AE - BE	0.06	0.04, 0.08	<0.001
BR - BE	-0.15	-0.19, -0.10	<0.001
BR - AE	-0.20	-0.26, -0.15	<0.001
AR - BE	-0.03	-0.04, -0.01	<0.001
AR - AE	-0.09	-0.11, -0.06	<0.001
AR - BR	0.12	0.07, 0.17	<0.001
Outcome
winner - loser	0.07	0.05, 0.09	<0.001
Time * Phase
Time * AE	-0.88	-1.1, -0.71	<0.001
Time * BR	0.02	-0.18, 0.21	0.9
Time * AR	-0.08	-0.11, -0.05	<0.001
Time * Outcome
Time * winner	-0.02	-0.03, -0.01	0.005
Phase * Outcome
AE * winner	0.05	0.03, 0.07	<0.001
BR * winner	0.02	-0.04, 0.08	0.5
AR * winner	0.14	0.12, 0.16	<0.001
Time * Phase * Outcome
Time * AE * winner	0.31	0.11, 0.50	0.002
Time * BR * winner	0.28	0.06, 0.50	0.013
Time * AR * winner	-0.04	-0.07, 0.00	0.033
Name.sd__(Intercept)	0.07
Name.cor__(Intercept).Time	-0.06
Name.sd__Time	0.06
Residual.sd__Observation	0.26
1 CI = Confidence Interval

Model Parameters

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

parameters::model_parameters(get(model))

# Fixed Effects

Parameter                              | Coefficient |       SE |         95% CI | t(169977) |      p
-----------------------------------------------------------------------------------------------------
(Intercept)                            |        0.52 | 4.45e-03 | [ 0.51,  0.52] |    115.70 | < .001
Time                                   |        0.06 | 5.60e-03 | [ 0.05,  0.07] |     10.19 | < .001
Phase [AE]                             |       -0.02 | 8.83e-03 | [-0.03,  0.00] |     -2.00 | 0.046 
Phase [BR]                             |       -0.15 |     0.03 | [-0.20, -0.09] |     -5.40 | < .001
Phase [AR]                             |       -0.10 | 9.83e-03 | [-0.12, -0.09] |    -10.65 | < .001
Outcome [winner]                       |        0.02 | 5.86e-03 | [ 0.01,  0.04] |      4.08 | < .001
Time × Phase [AE]                      |       -0.88 |     0.09 | [-1.05, -0.71] |    -10.21 | < .001
Time × Phase [BR]                      |        0.02 |     0.10 | [-0.18,  0.21] |      0.18 | 0.860 
Time × Phase [AR]                      |       -0.08 |     0.01 | [-0.11, -0.05] |     -5.24 | < .001
Time × Outcome [winner]                |       -0.02 | 7.20e-03 | [-0.03, -0.01] |     -2.83 | 0.005 
Phase [AE] × Outcome [winner]          |        0.05 |     0.01 | [ 0.03,  0.07] |      4.83 | < .001
Phase [BR] × Outcome [winner]          |        0.02 |     0.03 | [-0.04,  0.08] |      0.72 | 0.472 
Phase [AR] × Outcome [winner]          |        0.14 |     0.01 | [ 0.12,  0.16] |     12.52 | < .001
(Time × Phase [AE]) × Outcome [winner] |        0.31 |     0.10 | [ 0.11,  0.50] |      3.03 | 0.002 
(Time × Phase [BR]) × Outcome [winner] |        0.28 |     0.11 | [ 0.06,  0.50] |      2.49 | 0.013 
(Time × Phase [AR]) × Outcome [winner] |       -0.04 |     0.02 | [-0.07,  0.00] |     -2.14 | 0.033 

# Random Effects

Parameter                  | Coefficient
----------------------------------------
SD (Intercept: Name)       |        0.07
SD (Time: Name)            |        0.06
Cor (Intercept~Time: Name) |       -0.06
SD (Residual)              |        0.26

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.

FE: Time

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "0001"
terms=c("Time")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Phase

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "0002"
terms=c("Phase")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0


file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Phase

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "0003"
terms=c("Time", "Phase")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + lineR + line1 + stuff0 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Outcome

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "0004"
terms=c("Outcome")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Outcome

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "0005"
terms=c("Time", "Outcome")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + lineR + line1 + stuff0 + rect3

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Phase x Outcome

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "0006"
terms=c("Phase", "Outcome")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Outcome x Phase

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "0007"
terms=c("Time", "Outcome", "Phase")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 


gg88 <- gg88 + line0 + line1 + time0 + stuff0 + rect0

## gg88 <- gg88 + ggplot2::ylim(0.15,0.85)
## gg88 <- gg88 + ggplot2::coord_cartesian(ylim = c(0.25, 0.75))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

EF: Time:

NULL: Time

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "1001"
terms <- c("Time")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.54 | 0.52, 0.55
-0.50 |      0.53 | 0.52, 0.54
 0.00 |      0.52 | 0.52, 0.53
 0.50 |      0.52 | 0.51, 0.52
 1.00 |      0.51 | 0.50, 0.52

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Slope |       95% CI |     p
----------------------------
-0.02 | -0.03,  0.00 | 0.006

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05,)

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 
## gg88 <- gg88 + coord_cartesian(ylim = c(-0.25, 1.25))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

EF: Phase

NULL: Phase

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "1002"
terms <- c("Phase")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Phase | Predicted |     95% CI
------------------------------
BE    |      0.52 | 0.52, 0.53
AE    |      0.60 | 0.58, 0.62
BR    |      0.39 | 0.36, 0.43
AR    |      0.52 | 0.51, 0.53

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

Phase | Predicted |     95% CI |      p
---------------------------------------
BE    |      0.52 | 0.52, 0.53 | < .001
AE    |      0.60 | 0.58, 0.62 | < .001
BR    |      0.39 | 0.36, 0.43 | < .001
AR    |      0.52 | 0.51, 0.53 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Phase

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(
  limit_range=FALSE,  
  show_data = FALSE, 
  dot_alpha = 0.05)

gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Phase

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# Pairwise comparisons

Phase | Contrast |       95% CI |      p
----------------------------------------
BE-AE |    -0.08 | -0.09, -0.06 | < .001
BE-BR |     0.13 |  0.10,  0.17 | < .001
BE-AR | 4.18e-03 | -0.01,  0.01 | 0.442 
AE-BR |     0.21 |  0.17,  0.25 | < .001
AE-AR |     0.08 |  0.06,  0.10 | < .001
BR-AR |    -0.13 | -0.17, -0.09 | < .001

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Phase

NULL: Time x Phase

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "1003"
terms <- c("Phase", "Time")
terms <- c("Time", "Phase")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Phase: BE

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.48 | 0.48,  0.49
-0.50 |      0.51 | 0.50,  0.51
 0.00 |      0.53 | 0.52,  0.53
 0.50 |      0.55 | 0.54,  0.56
 1.00 |      0.57 | 0.56,  0.58

Phase: AE

 Time | Predicted |       95% CI
--------------------------------
-1.00 |      1.17 |  1.08,  1.27
-0.50 |      0.86 |  0.81,  0.91
 0.00 |      0.55 |  0.54,  0.56
 0.50 |      0.23 |  0.19,  0.27
 1.00 |     -0.08 | -0.16,  0.00

Phase: BR

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.14 | 0.02,  0.27
-0.50 |      0.27 | 0.20,  0.35
 0.00 |      0.40 | 0.37,  0.42
 0.50 |      0.53 | 0.50,  0.55
 1.00 |      0.65 | 0.58,  0.72

Phase: AR

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.57 | 0.55,  0.60
-0.50 |      0.55 | 0.53,  0.56
 0.00 |      0.52 | 0.51,  0.53
 0.50 |      0.49 | 0.49,  0.50
 1.00 |      0.46 | 0.46,  0.47

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Phase | Slope |       95% CI |      p
-------------------------------------
BE    |  0.06 |  0.05,  0.07 | < .001
AE    | -0.83 | -0.99, -0.66 | < .001
BR    |  0.07 | -0.12,  0.27 | 0.452 
AR    | -0.02 | -0.05,  0.01 | 0.219 

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time x Phase

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
ggsave(file = "file.svg", plot = gg88, width=8, height=6)
gg88

DIFF: Time x Phase

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# (Average) Linear trend for Time

Phase | Contrast |       95% CI |      p
----------------------------------------
BE-AE |     0.88 |  0.71,  1.05 | < .001
BE-BR |    -0.02 | -0.21,  0.18 | 0.860 
BE-AR |     0.08 |  0.05,  0.11 | < .001
AE-BR |    -0.90 | -1.16, -0.64 | < .001
AE-AR |    -0.81 | -0.98, -0.63 | < .001
BR-AR |     0.09 | -0.10,  0.29 | 0.414 

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Outcome

NULL: Outcome

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "1004"
terms <- c("Outcome")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Outcome | Predicted |     95% CI
--------------------------------
loser   |      0.45 | 0.44, 0.46
winner  |      0.52 | 0.51, 0.52

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

Outcome | Predicted |     95% CI |      p
-----------------------------------------
loser   |      0.45 | 0.44, 0.46 | < .001
winner  |      0.52 | 0.51, 0.52 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Outcome

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
ggsave(file = "file.svg", plot = gg88, width=8, height=6)
gg88

DIFF: Outcome

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# Pairwise comparisons

Outcome      | Contrast |       95% CI |      p
-----------------------------------------------
loser-winner |    -0.07 | -0.08, -0.06 | < .001

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Outcome:

NULL: Time x Outcome

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "1005"
terms <- c("Time", "Outcome")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Outcome: loser

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.50 | 0.47, 0.52
-0.50 |      0.48 | 0.47, 0.50
 0.00 |      0.47 | 0.46, 0.48
 0.50 |      0.46 | 0.45, 0.47
 1.00 |      0.45 | 0.43, 0.47

Outcome: winner

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.55 | 0.54, 0.57
-0.50 |      0.55 | 0.54, 0.56
 0.00 |      0.54 | 0.53, 0.55
 0.50 |      0.53 | 0.53, 0.54
 1.00 |      0.53 | 0.51, 0.54

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Outcome | Slope |     95% CI |      p
-------------------------------------
loser   |  0.06 | 0.05, 0.07 | < .001
winner  |  0.04 | 0.03, 0.05 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT Time x Outcome

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + line0 + line1 + time0 + stuff0 + rect3

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Time x Outcome

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# (Average) Linear trend for Time

Outcome      | Contrast |     95% CI |     p
--------------------------------------------
loser-winner |     0.02 | 0.01, 0.03 | 0.005

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Phase x Outcome

NULL: Phase x Outcome

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "1006"
terms <- c("Phase", "Outcome")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Outcome: loser

Phase | Predicted |     95% CI
------------------------------
BE    |      0.51 | 0.50, 0.52
AE    |      0.55 | 0.52, 0.58
BR    |      0.36 | 0.30, 0.43
AR    |      0.41 | 0.39, 0.43

Outcome: winner

Phase | Predicted |     95% CI
------------------------------
BE    |      0.53 | 0.53, 0.54
AE    |      0.61 | 0.59, 0.62
BR    |      0.39 | 0.35, 0.43
AR    |      0.58 | 0.57, 0.59

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

Phase | Outcome | Predicted |     95% CI |      p
-------------------------------------------------
BE    |   loser |      0.51 | 0.50, 0.52 | < .001
AE    |   loser |      0.55 | 0.52, 0.58 | < .001
BR    |   loser |      0.36 | 0.30, 0.43 | < .001
AR    |   loser |      0.41 | 0.39, 0.43 | < .001
BE    |  winner |      0.53 | 0.53, 0.54 | < .001
AE    |  winner |      0.61 | 0.59, 0.62 | < .001
BR    |  winner |      0.39 | 0.35, 0.43 | < .001
AR    |  winner |      0.58 | 0.57, 0.59 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Phase x Outcome

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)

gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

PLOT: Phase x Outcome [low-level]

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "1006"
terms <- c("Phase", "Outcome")

type <- "eff"
type <- "pred"

suppressWarnings(rm(list = ls(pattern = "^gg44")))
gg44 <- sjPlot::plot_model(
  get(model), 
  type = type, 
  terms = terms, 
  )

gg44 <- gg44 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjplot-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg44, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjplot-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg44, width=8, height=6)
gg44

DIFF: Phase x Outcome

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# Pairwise comparisons

Phase |       Outcome | Contrast |       95% CI |      p
--------------------------------------------------------
BE-AE |   loser-loser |    -0.04 | -0.07, -0.02 | 0.003 
BE-BR |   loser-loser |     0.15 |  0.08,  0.21 | < .001
BE-AR |   loser-loser |     0.10 |  0.08,  0.12 | < .001
BE-BE |  loser-winner |    -0.03 | -0.04, -0.01 | < .001
BE-AE |  loser-winner |    -0.10 | -0.12, -0.08 | < .001
BE-BR |  loser-winner |     0.12 |  0.08,  0.16 | < .001
BE-AR |  loser-winner |    -0.07 | -0.08, -0.05 | < .001
AE-BR |   loser-loser |     0.19 |  0.12,  0.26 | < .001
AE-AR |   loser-loser |     0.14 |  0.11,  0.18 | < .001
AE-BE |  loser-winner |     0.02 | -0.01,  0.05 | 0.265 
AE-AE |  loser-winner |    -0.06 | -0.09, -0.02 | 0.001 
AE-BR |  loser-winner |     0.16 |  0.11,  0.21 | < .001
AE-AR |  loser-winner |    -0.03 | -0.06,  0.00 | 0.103 
BR-AR |   loser-loser |    -0.05 | -0.12,  0.02 | 0.201 
BR-BE |  loser-winner |    -0.17 | -0.24, -0.11 | < .001
BR-AE |  loser-winner |    -0.24 | -0.31, -0.18 | < .001
BR-BR |  loser-winner |    -0.03 | -0.10,  0.05 | 0.463 
BR-AR |  loser-winner |    -0.22 | -0.28, -0.15 | < .001
AR-BE |  loser-winner |    -0.12 | -0.15, -0.10 | < .001
AR-AE |  loser-winner |    -0.20 | -0.23, -0.17 | < .001
AR-BR |  loser-winner |     0.02 | -0.02,  0.06 | 0.400 
AR-AR |  loser-winner |    -0.17 | -0.19, -0.14 | < .001
BE-AE | winner-winner |    -0.07 | -0.09, -0.06 | < .001
BE-BR | winner-winner |     0.14 |  0.11,  0.18 | < .001
BE-AR | winner-winner |    -0.04 | -0.06, -0.03 | < .001
AE-BR | winner-winner |     0.22 |  0.18,  0.26 | < .001
AE-AR | winner-winner |     0.03 |  0.01,  0.05 | 0.006 
BR-AR | winner-winner |    -0.19 | -0.22, -0.15 | < .001

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Phase x Outcome

NULL: Time x Phase x Outcome

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

extra <- "1007"
terms <- c("Time", "Phase", "Outcome")
terms <- c("Time", "Outcome", "Phase") ## CAUTION
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Outcome: loser
Phase: BE

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.45 | 0.44,  0.46
-0.50 |      0.48 | 0.47,  0.49
 0.00 |      0.51 | 0.50,  0.52
 0.50 |      0.54 | 0.53,  0.55
 1.00 |      0.57 | 0.56,  0.59

Outcome: loser
Phase: AE

 Time | Predicted |       95% CI
--------------------------------
-1.00 |      1.32 |  1.13,  1.50
-0.50 |      0.91 |  0.81,  1.01
 0.00 |      0.49 |  0.48,  0.51
 0.50 |      0.08 |  0.01,  0.16
 1.00 |     -0.33 | -0.48, -0.17

Outcome: loser
Phase: BR

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.29 | 0.04,  0.54
-0.50 |      0.33 | 0.18,  0.48
 0.00 |      0.37 | 0.31,  0.42
 0.50 |      0.41 | 0.36,  0.45
 1.00 |      0.44 | 0.30,  0.59

Outcome: loser
Phase: AR

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.43 | 0.38,  0.47
-0.50 |      0.42 | 0.38,  0.45
 0.00 |      0.41 | 0.39,  0.43
 0.50 |      0.40 | 0.39,  0.41
 1.00 |      0.39 | 0.37,  0.41

Outcome: winner
Phase: BE

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.50 | 0.49,  0.51
-0.50 |      0.52 | 0.51,  0.52
 0.00 |      0.54 | 0.53,  0.54
 0.50 |      0.56 | 0.55,  0.57
 1.00 |      0.58 | 0.56,  0.59

Outcome: winner
Phase: AE

 Time | Predicted |       95% CI
--------------------------------
-1.00 |      1.11 |  0.99,  1.22
-0.50 |      0.84 |  0.78,  0.90
 0.00 |      0.57 |  0.56,  0.58
 0.50 |      0.30 |  0.26,  0.34
 1.00 |      0.03 | -0.06,  0.13

Outcome: winner
Phase: BR

 Time | Predicted |       95% CI
--------------------------------
-1.00 |      0.08 | -0.06,  0.21
-0.50 |      0.24 |  0.16,  0.33
 0.00 |      0.41 |  0.38,  0.44
 0.50 |      0.58 |  0.56,  0.61
 1.00 |      0.75 |  0.67,  0.83

Outcome: winner
Phase: AR

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.65 | 0.62,  0.67
-0.50 |      0.61 | 0.59,  0.63
 0.00 |      0.57 | 0.56,  0.58
 0.50 |      0.54 | 0.53,  0.54
 1.00 |      0.50 | 0.49,  0.51

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Outcome | Phase | Slope |       95% CI |      p
-----------------------------------------------
loser   |    BE |  0.06 |  0.05,  0.07 | < .001
loser   |    AE | -0.83 | -0.99, -0.66 | < .001
loser   |    BR |  0.07 | -0.12,  0.27 | 0.452 
loser   |    AR | -0.02 | -0.05,  0.01 | 0.187 
winner  |    BE |  0.04 |  0.03,  0.05 | < .001
winner  |    AE | -0.54 | -0.64, -0.44 | < .001
winner  |    BR |  0.34 |  0.23,  0.44 | < .001
winner  |    AR | -0.08 | -0.09, -0.06 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time x Phase x Outcome

extra <- "1007"
terms <- c("Time", "Phase", "Outcome")
terms <- c("Time", "Outcome", "Phase") ## CAUTION
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 

gg88 <- gg88 + rect0

# gg88 <- gg88 + ggplot2::ylim(0.10,0.90)
gg88 <- gg88 + scale_y_continuous(breaks = seq(0,1,by=0.1))

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=16)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=16)
gg88

knitr::opts_chunk$set(fig.width=unit(12,"cm"), fig.height=unit(18,"cm"))

gg88

PLOT: Time x Phase x Outcome (low-level)

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
          get(model),
          type="emm",
          terms=terms)

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 + rect0
## gg88 <- gg88 + ggplot2::ylim(0.15,0.85)

file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjPlot-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=4, height=16)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjPlot-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=4, height=16)
gg88

DIFF: Time x Phase x Outcome

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# (Average) Linear trend for Time

Outcome       | Phase | Contrast |       95% CI |      p
--------------------------------------------------------
loser-loser   | BE-AE |     0.88 |  0.71,  1.05 | < .001
loser-loser   | BE-BR |    -0.02 | -0.21,  0.18 | 0.860 
loser-loser   | BE-AR |     0.08 |  0.05,  0.11 | < .001
loser-winner  | BE-BE |     0.02 |  0.01,  0.03 | 0.006 
loser-winner  | BE-AE |     0.60 |  0.49,  0.70 | < .001
loser-winner  | BE-BR |    -0.28 | -0.39, -0.17 | < .001
loser-winner  | BE-AR |     0.13 |  0.11,  0.15 | < .001
loser-loser   | AE-BR |    -0.90 | -1.16, -0.64 | < .001
loser-loser   | AE-AR |    -0.81 | -0.98, -0.63 | < .001
loser-winner  | AE-BE |    -0.86 | -1.03, -0.69 | < .001
loser-winner  | AE-AE |    -0.29 | -0.48, -0.09 | 0.006 
loser-winner  | AE-BR |    -1.16 | -1.36, -0.96 | < .001
loser-winner  | AE-AR |    -0.75 | -0.92, -0.58 | < .001
loser-loser   | BR-AR |     0.09 | -0.10,  0.29 | 0.371 
loser-winner  | BR-BE |     0.04 | -0.16,  0.23 | 0.729 
loser-winner  | BR-AE |     0.61 |  0.39,  0.83 | < .001
loser-winner  | BR-BR |    -0.26 | -0.48, -0.04 | 0.024 
loser-winner  | BR-AR |     0.15 | -0.04,  0.35 | 0.145 
loser-winner  | AR-BE |    -0.06 | -0.09, -0.03 | < .001
loser-winner  | AR-AE |     0.52 |  0.41,  0.63 | < .001
loser-winner  | AR-BR |    -0.36 | -0.47, -0.25 | < .001
loser-winner  | AR-AR |     0.06 |  0.02,  0.09 | < .001
winner-winner | BE-AE |     0.58 |  0.47,  0.68 | < .001
winner-winner | BE-BR |    -0.30 | -0.41, -0.19 | < .001
winner-winner | BE-AR |     0.11 |  0.10,  0.13 | < .001
winner-winner | AE-BR |    -0.87 | -1.02, -0.73 | < .001
winner-winner | AE-AR |    -0.46 | -0.57, -0.36 | < .001
winner-winner | BR-AR |     0.41 |  0.30,  0.52 | < .001

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

CHECK: Model Assumptions

pp(model);

fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

performance::check_model(get(model))

ggsave(
  file = "performance-test.png", 
  ## file.path(dir8, "summary-combined-models.perf0.compare-02-score.png"),
  ## plot = gg88,
  ## plot = last_plot(),
  width=8,
  height=12)

Warning: Failed to fit group -1.
Failed to fit group -1.
Caused by error in `predLoess()`:
! workspace required (43350042532) is too large probably because of setting 'se = TRUE'.

ggsave(
  file = "performance-test.png", 
  ## file.path(dir8, "summary-combined-models.perf0.compare-02-score.png"),
  ## plot = gg88,
  ## plot = last_plot(),
  width=8,
  height=12)

Warning: Failed to fit group -1.
Failed to fit group -1.
Caused by error in `predLoess()`:
! workspace required (43350042532) is too large probably because of setting 'se = TRUE'.

ls()

 [1] "cat0"                     "check_pkgs"               "control"                  "create_factor_mappings"   "df0"                      "df2"                      "df3"                      "dir8"                    
 [9] "extra"                    "fbase"                    "file"                     "fit01aPh"                 "fit02aPh"                 "fit03aPh"                 "fit04aPh"                 "fpath"                   
[17] "generate_permutations"    "get_stars"                "gg44"                     "gg88"                     "ggpairs_lower_fun"        "hasConverged"             "ifd0"                     "ifn0"                    
[25] "jiko_count_na"            "line0"                    "line1"                    "line2"                    "lineR"                    "model"                    "name"                     "obj_has_name"            
[33] "ofd0"                     "ofd2"                     "pp"                       "pred0"                    "pred0tab"                 "rect0"                    "rect2"                    "rect3"                   
[41] "REML"                     "replace_factor_levels"    "replace_outliers_with_NA" "rescale01"                "rxy0"                     "rxy2"                     "rxy3"                     "sep0"                    
[49] "sep1"                     "sep2"                     "set.w"                    "stuff0"                   "tbl0"                     "terms"                    "test0"                    "test0tab"                
[57] "test2"                    "test2tab"                 "time0"                    "type"                     "var0"                     "write_model_info"        

Auxiliary Checkups

Variables influence / importance in fit04aPh

## TODO CHANGE TITLE
formula(fit04aPh)

Agency ~ (Time | Name) + Time * Phase * Outcome

eff_size0 <- effectsize::eta_squared(fit04aPh)
eff_size2 <- eff_size0 %>%
  dplyr::mutate(Interpret = effectsize::interpret_eta_squared(Eta2_partial))

performance::print_html(eff_size2)

Effect Size for ANOVA
Parameter	Eta2 (partial)	95% CI	Interpret
Time	2.41e-04	[0.00, 1.00]	very small
Phase	6.86e-04	[0.00, 1.00]	very small
Outcome	5.68e-03	[0.00, 1.00]	very small
Time:Phase	2.04e-03	[0.00, 1.00]	very small
Time:Outcome	5.84e-05	[0.00, 1.00]	very small
Phase:Outcome	1.01e-03	[0.00, 1.00]	very small
Time:Phase:Outcome	1.24e-04	[0.00, 1.00]	very small

List models (for Model Hierarchy / Hierarchies)

temp <- lapply(ls(pattern="^fit\\d+"), pp)

fit01aPh: [df0] Agency ~ (1 | Name) + 1
fit02aPh: [df0] Agency ~ (Time | Name) + Time
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

Compare models using performance package with score

perf0 <- performance::compare_performance(
  fit01aPh, # [df0] Agency ~ (1 | Name) + 1
  fit02aPh, # [df0] Agency ~ (Time | Name) + Time
  fit03aPh, # [df0] Agency ~ (Time | Name) + Time * Phase
  fit04aPh, # [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
  ## CAUTION: COMMA
  rank = TRUE, verbose = FALSE)
perf0 %>% performance::print_html()

Comparison of Model Performance Indices
Name	Model	R2 (cond.)	R2 (marg.)	ICC	RMSE	Sigma	AIC weights	AICc weights	BIC weights	Performance-Score
fit04aPh	lmerModLmerTest	0.10	0.02	0.08	0.26	0.26	1.00	1.00	1.00	84.49%
fit03aPh	lmerModLmerTest	0.11	5.28e-03	0.10	0.26	0.26	4.39e-118	4.39e-118	1.23e-100	51.83%
fit02aPh	lmerModLmerTest	0.10	4.62e-04	0.10	0.26	0.26	2.04e-289	2.04e-289	6.95e-259	41.69%
fit01aPh	lmerModLmerTest	0.08	0.00	0.08	0.26	0.26	0.00e+00	0.00e+00	0.00e+00	1.00%
NA

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- perf0 %>% plot()
ggsave(
  file = file.path(dir8, "summary-combined-models.perf0.compare-02-score.png"),
  plot = gg88,
  width=8,
  height=6)

gg88

Compare less models

perf2 <- performance::compare_performance(
  ## fit01aPh, # [df0] Agency ~ (1 | Name) + 1
  ## fit02aPh, # [df0] Agency ~ (Time | Name) + Time
  fit03aPh, # [df0] Agency ~ (Time | Name) + Time * Phase
  fit04aPh, # [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
  ## CAUTION: COMMA
  rank = TRUE, verbose = FALSE)
perf2 %>% performance::print_html()

Comparison of Model Performance Indices
Name	Model	R2 (cond.)	R2 (marg.)	ICC	RMSE	Sigma	AIC weights	AICc weights	BIC weights	Performance-Score
fit04aPh	lmerModLmerTest	0.10	0.02	0.08	0.26	0.26	1.00	1.00	1.00	75.00%
fit03aPh	lmerModLmerTest	0.11	5.28e-03	0.10	0.26	0.26	4.39e-118	4.39e-118	1.23e-100	25.00%
NA

suppressWarnings(rm(list = ls(pattern = "^gg44")))
gg44 <- perf2 %>% plot()
ggsave(
  file = file.path(dir8, "summary-combined-models.perf2.compare-03-score.png"),
  plot = gg44,
  width=8,
  height=6)

gg44  

Comment

Add interpretation of the difference between best fit to data and best predictive power and overfitting (aspecially in the context of the last model).

gg44

Performance Table Sorted by R2_conditional

perf0 %>% dplyr::arrange(desc(R2_conditional)) %>% performance::print_html()

Comparison of Model Performance Indices
Name	Model	R2 (cond.)	R2 (marg.)	ICC	RMSE	Sigma	AIC weights	AICc weights	BIC weights	Performance-Score
fit03aPh	lmerModLmerTest	0.11	5.28e-03	0.10	0.26	0.26	4.39e-118	4.39e-118	1.23e-100	51.83%
fit02aPh	lmerModLmerTest	0.10	4.62e-04	0.10	0.26	0.26	2.04e-289	2.04e-289	6.95e-259	41.69%
fit04aPh	lmerModLmerTest	0.10	0.02	0.08	0.26	0.26	1.00	1.00	1.00	84.49%
fit01aPh	lmerModLmerTest	0.08	0.00	0.08	0.26	0.26	0.00e+00	0.00e+00	0.00e+00	1.00%
NA

Performance Table Sorted by R2_marginal

perf0 %>% dplyr::arrange(desc(R2_marginal)) %>% performance::print_html()

Comparison of Model Performance Indices
Name	Model	R2 (cond.)	R2 (marg.)	ICC	RMSE	Sigma	AIC weights	AICc weights	BIC weights	Performance-Score
fit04aPh	lmerModLmerTest	0.10	0.02	0.08	0.26	0.26	1.00	1.00	1.00	84.49%
fit03aPh	lmerModLmerTest	0.11	5.28e-03	0.10	0.26	0.26	4.39e-118	4.39e-118	1.23e-100	51.83%
fit02aPh	lmerModLmerTest	0.10	4.62e-04	0.10	0.26	0.26	2.04e-289	2.04e-289	6.95e-259	41.69%
fit01aPh	lmerModLmerTest	0.08	0.00	0.08	0.26	0.26	0.00e+00	0.00e+00	0.00e+00	1.00%
NA

Interpret R2

model <- "fit01aPh" # [df0] Agency ~ (1 | Name) + 1
model <- "fit02aPh" # [df0] Agency ~ (Time | Name) + Time
model <- "fit03aPh" # [df0] Agency ~ (Time | Name) + Time * Phase
model <- "fit04aPh" # [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

cat0(effectsize::interpret_r2(performance::r2(get(model))$R2_conditional, rules="cohen1988"))

weak 

cat0(effectsize::interpret_r2(performance::r2(get(model))$R2_marginal, rules="cohen1988"))

weak 

Plot models: Part 1: Basic Models

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_models(
  ## CAUTION the null model can not be used here 
  ## Thus to keep the numbers consistent I have 
  ## used model 02 as an input twice
  ## fit01aPh, # [df0] Agency ~ (Time | Name) + Time
  fit02aPh, # [df0] Agency ~ (Time | Name) + Time
  fit03aPh, # [df0] Agency ~ (Time | Name) + Time * Phase
  fit04aPh, # [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
  spacing=1, dot.size=1
)

ggsave(
  file = file.path(dir8, "summary-all-models-part-01.png"),
  plot = gg88,
  width=5,
  height=4)

gg88

Tab Models

The above `plot_models` seems to be a lot more illustrative

# , eval = TRUE, results='hide'
library(sjPlot)

Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!

library(sjmisc)

Attaching package: 'sjmisc'

The following object is masked from 'package:skimr':

    to_long

The following object is masked from 'package:purrr':

    is_empty

The following object is masked from 'package:tidyr':

    replace_na

The following object is masked from 'package:tibble':

    add_case

library(sjlabelled)

Attaching package: 'sjlabelled'

The following object is masked from 'package:forcats':

    as_factor

The following object is masked from 'package:dplyr':

    as_label

The following object is masked from 'package:ggplot2':

    as_label

tab0 <- sjPlot::tab_model(
  fit01aPh, # [df0] Agency ~ (1 | Name) + 1
  fit02aPh, # [df0] Agency ~ (Time | Name) + Time
  fit03aPh, # [df0] Agency ~ (Time | Name) + Time * Phase
  fit04aPh, # [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
    show.reflvl = TRUE,
    show.intercept = TRUE,
    p.style = "numeric_stars",
    file = file.path(dir8, "summary-all-models-part-01-SEE-PNG-images-for-a-better-overview.html"))
tab0

	Agency			Agency			Agency			Agency
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	0.50 ***	0.49 – 0.50	<0.001	0.50 ***	0.49 – 0.50	<0.001	0.53 ***	0.52 – 0.53	<0.001	0.52 ***	0.51 – 0.52	<0.001
BE	Reference			Reference			Reference			Reference
AE							0.02 **	0.01 – 0.02	0.002	-0.02 *	-0.03 – -0.00	0.046
BR							-0.13 ***	-0.16 – -0.11	<0.001	-0.15 ***	-0.20 – -0.09	<0.001
PhaseAE:Outcomewinner										0.05 ***	0.03 – 0.07	<0.001
AR							-0.00	-0.01 – 0.01	0.899	-0.10 ***	-0.12 – -0.09	<0.001
PhaseAR:Outcomewinner										0.14 ***	0.12 – 0.16	<0.001
loser	Reference			Reference			Reference			Reference
winner										0.02 ***	0.01 – 0.04	<0.001
PhaseBR:Outcomewinner										0.02	-0.04 – 0.08	0.472
Time				-0.01 ***	-0.02 – -0.01	<0.001	0.04 ***	0.04 – 0.05	<0.001	0.06 ***	0.05 – 0.07	<0.001
Time:Outcomewinner										-0.02 **	-0.03 – -0.01	0.005
Time:PhaseAE							-0.64 ***	-0.73 – -0.55	<0.001	-0.88 ***	-1.05 – -0.71	<0.001
Time:PhaseAE:Outcomewinner										0.31 **	0.11 – 0.50	0.002
Time:PhaseAR							-0.11 ***	-0.12 – -0.10	<0.001	-0.08 ***	-0.11 – -0.05	<0.001
Time:PhaseAR:Outcomewinner										-0.04 *	-0.07 – -0.00	0.033
Time:PhaseBR							0.22 ***	0.13 – 0.32	<0.001	0.02	-0.18 – 0.21	0.860
Time:PhaseBR:Outcomewinner										0.28 *	0.06 – 0.50	0.013
Random Effects
σ2	0.07			0.07			0.07			0.07
τ00	0.01 Name			0.01 Name			0.01 Name			0.00 Name
τ11				0.00 Name.Time			0.00 Name.Time			0.00 Name.Time
ρ01				0.18 Name			0.18 Name			-0.06 Name
ICC	0.08			0.10			0.10			0.08
N	870 Name			870 Name			870 Name			870 Name
Observations	169997			169997			169997			169997
Marginal R2 / Conditional R2	0.000 / 0.084			0.000 / 0.102			0.005 / 0.106			0.020 / 0.101
* p<0.05   ** p<0.01   *** p<0.001

## knitr::html
## papaja::
## rmarkdown::render
## rmarkdown::render(tab0$knitr)
## tab0$knitr

Tabulate models (publication ready)

tab0 <- sjPlot::tab_model(
  fit01aPh, # [df0] Agency ~ (1 | Name) + 1
  fit02aPh, # [df0] Agency ~ (Time | Name) + Time
  fit03aPh, # [df0] Agency ~ (Time | Name) + Time * Phase
  fit04aPh, # [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
    show.reflvl = FALSE,
    show.intercept = TRUE,
    p.style = "numeric_stars",
    file = file.path(dir8, "summary-all-models-part-02-publication.html"))
tab0

	Agency			Agency			Agency			Agency
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	0.50 ***	0.49 – 0.50	<0.001	0.50 ***	0.49 – 0.50	<0.001	0.53 ***	0.52 – 0.53	<0.001	0.52 ***	0.51 – 0.52	<0.001
Time				-0.01 ***	-0.02 – -0.01	<0.001	0.04 ***	0.04 – 0.05	<0.001	0.06 ***	0.05 – 0.07	<0.001
Phase [AE]							0.02 **	0.01 – 0.02	0.002	-0.02 *	-0.03 – -0.00	0.046
Phase [BR]							-0.13 ***	-0.16 – -0.11	<0.001	-0.15 ***	-0.20 – -0.09	<0.001
Phase [AR]							-0.00	-0.01 – 0.01	0.899	-0.10 ***	-0.12 – -0.09	<0.001
Time × Phase [AE]							-0.64 ***	-0.73 – -0.55	<0.001	-0.88 ***	-1.05 – -0.71	<0.001
Time × Phase [BR]							0.22 ***	0.13 – 0.32	<0.001	0.02	-0.18 – 0.21	0.860
Time × Phase [AR]							-0.11 ***	-0.12 – -0.10	<0.001	-0.08 ***	-0.11 – -0.05	<0.001
Outcome [winner]										0.02 ***	0.01 – 0.04	<0.001
Time × Outcome [winner]										-0.02 **	-0.03 – -0.01	0.005
Phase [AE] × Outcome [winner]										0.05 ***	0.03 – 0.07	<0.001
Phase [BR] × Outcome [winner]										0.02	-0.04 – 0.08	0.472
Phase [AR] × Outcome [winner]										0.14 ***	0.12 – 0.16	<0.001
(Time × Phase [AE]) × Outcome [winner]										0.31 **	0.11 – 0.50	0.002
(Time × Phase [BR]) × Outcome [winner]										0.28 *	0.06 – 0.50	0.013
(Time × Phase [AR]) × Outcome [winner]										-0.04 *	-0.07 – -0.00	0.033
Random Effects
σ2	0.07			0.07			0.07			0.07
τ00	0.01 Name			0.01 Name			0.01 Name			0.00 Name
τ11				0.00 Name.Time			0.00 Name.Time			0.00 Name.Time
ρ01				0.18 Name			0.18 Name			-0.06 Name
ICC	0.08			0.10			0.10			0.08
N	870 Name			870 Name			870 Name			870 Name
Observations	169997			169997			169997			169997
Marginal R2 / Conditional R2	0.000 / 0.084			0.000 / 0.102			0.005 / 0.106			0.020 / 0.101
* p<0.05   ** p<0.01   *** p<0.001

# save.image(file=file.path(dir8, "session.RData"))
# load(file.path(dir8, "session.RData")

Report

result <- report::report(fit04aPh)

print(result)

We fitted a linear mixed model (estimated using REML and Nelder-Mead optimizer) to predict Agency with Time, Phase and Outcome (formula: Agency ~ Time * Phase * Outcome). The model included Time as random effects
(formula: ~Time | Name). The model's total explanatory power is weak (conditional R2 = 0.10) and the part related to the fixed effects alone (marginal R2) is of 0.02. The model's intercept, corresponding to Time = 0,
Phase = BE and Outcome = loser, is at 0.52 (95% CI [0.51, 0.52], t(169977) = 115.70, p < .001). Within this model:

  - The effect of Time is statistically significant and positive (beta = 0.06, 95% CI [0.05, 0.07], t(169977) = 10.19, p < .001; Std. beta = 0.12, 95% CI [0.10, 0.14])
  - The effect of Phase [AE] is statistically significant and negative (beta = -0.02, 95% CI [-0.03, -3.46e-04], t(169977) = -2.00, p = 0.046; Std. beta = 0.16, 95% CI [0.06, 0.26])
  - The effect of Phase [BR] is statistically significant and negative (beta = -0.15, 95% CI [-0.20, -0.09], t(169977) = -5.40, p < .001; Std. beta = -0.54, 95% CI [-0.78, -0.30])
  - The effect of Phase [AR] is statistically significant and negative (beta = -0.10, 95% CI [-0.12, -0.09], t(169977) = -10.65, p < .001; Std. beta = -0.37, 95% CI [-0.44, -0.29])
  - The effect of Outcome [winner] is statistically significant and positive (beta = 0.02, 95% CI [0.01, 0.04], t(169977) = 4.08, p < .001; Std. beta = 0.09, 95% CI [0.05, 0.13])
  - The effect of Time × Phase [AE] is statistically significant and negative (beta = -0.88, 95% CI [-1.05, -0.71], t(169977) = -10.21, p < .001; Std. beta = -1.86, 95% CI [-2.22, -1.51])
  - The effect of Time × Phase [BR] is statistically non-significant and positive (beta = 0.02, 95% CI [-0.18, 0.21], t(169977) = 0.18, p = 0.860; Std. beta = 0.04, 95% CI [-0.37, 0.45])
  - The effect of Time × Phase [AR] is statistically significant and negative (beta = -0.08, 95% CI [-0.11, -0.05], t(169977) = -5.24, p < .001; Std. beta = -0.16, 95% CI [-0.22, -0.10])
  - The effect of Time × Outcome [winner] is statistically significant and negative (beta = -0.02, 95% CI [-0.03, -6.25e-03], t(169977) = -2.83, p = 0.005; Std. beta = -0.04, 95% CI [-0.07, -0.01])
  - The effect of Phase [AE] × Outcome [winner] is statistically significant and positive (beta = 0.05, 95% CI [0.03, 0.07], t(169977) = 4.83, p < .001; Std. beta = 0.11, 95% CI [-8.07e-03, 0.23])
  - The effect of Phase [BR] × Outcome [winner] is statistically non-significant and positive (beta = 0.02, 95% CI [-0.04, 0.08], t(169977) = 0.72, p = 0.472; Std. beta = 0.01, 95% CI [-0.27, 0.29])
  - The effect of Phase [AR] × Outcome [winner] is statistically significant and positive (beta = 0.14, 95% CI [0.12, 0.16], t(169977) = 12.52, p < .001; Std. beta = 0.53, 95% CI [0.44, 0.62])
  - The effect of (Time × Phase [AE]) × Outcome [winner] is statistically significant and positive (beta = 0.31, 95% CI [0.11, 0.50], t(169977) = 3.03, p = 0.002; Std. beta = 0.65, 95% CI [0.23, 1.07])
  - The effect of (Time × Phase [BR]) × Outcome [winner] is statistically significant and positive (beta = 0.28, 95% CI [0.06, 0.50], t(169977) = 2.49, p = 0.013; Std. beta = 0.59, 95% CI [0.13, 1.06])
  - The effect of (Time × Phase [AR]) × Outcome [winner] is statistically significant and negative (beta = -0.04, 95% CI [-0.07, -2.93e-03], t(169977) = -2.14, p = 0.033; Std. beta = -0.08, 95% CI [-0.14, -6.19e-03])

Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using a Wald t-distribution approximation.